package biolabs.model.react

import scala.collection.mutable.{Map => MMap}
import scala.math.{max, pow}

object Reaction {
	def react(
			concentrations: Map[Any, Double],
			reactions: Iterable[Reaction]
	): Map[Any, Double] = {
		val oneWays = reactions.foldLeft(List[OneWayReaction]())(
			(list, reaction) =>
				if(reaction.backward != null)
					reaction.forward :: reaction.backward :: list
				else
					reaction.forward :: list
		)
		return performReactions(concentrations, oneWays);
	}

	private def performReactions(
			concentrations: Map[Any, Double],
			reactions: Iterable[OneWayReaction]
	): Map[Any, Double] = {
		val actualRates = computeActualRates(concentrations, reactions)
		val conc = MMap() ++ concentrations
		for((reaction, actualRate) <- actualRates) {
			for((factor, sub) <- reaction.factors)
				conc(sub) = max(0, conc(sub) + factor * actualRate)
		}
		return conc.toMap
	}

	private def computeActualRates(
			concentrations: Map[Any, Double],
			reactions: Iterable[OneWayReaction]
	): Iterable[(OneWayReaction, Double)] = {
		val desiredRates = computeDesiredRates(concentrations, reactions)
		return _computeActualRates(concentrations, desiredRates)
	}

	private def _computeActualRates(
			concentrations: Map[Any, Double],
			desiredRates: Iterable[(OneWayReaction, Double)]
	): Iterable[(OneWayReaction, Double)] = {

		// compute concentrations that would satisfy all reactions
		val desiredConc = MMap() ++ concentrations.map { case (sub, _) => (sub, 0.0) }
		for((reaction, desiredRate) <- desiredRates) {
			for((factor, sub) <- reaction.reactants)
				desiredConc(sub) += factor * desiredRate
		}

		// identify the most limiting substance minS
		// and its satisfaction ratio minQ
		var minS: Any = null
		var minQ = Double.PositiveInfinity
		for((sub, c) <- concentrations) {
			val q = c / desiredConc(sub)
			if(q <= minQ) {
				minS = sub
				minQ = q
			}
		}

		if(minQ >= 1)
			return desiredRates

		// fulfill all reactions at ratio minQ
		val actualRates = desiredRates.map { case (reaction, rate) => (reaction, rate * minQ) }

		// compute the remaining concentrations, exclude minS and all concentrations <= 0
		var excluded = List(minS)
		val remainingConc = (concentrations.map { case (sub, c) => (sub, c - minQ*desiredConc(sub)) } - minS).filter {
			case (sub, c) => if(c > 0) true else {excluded = sub :: excluded; false}
		}

		// recompute desired rates, exclude all reactions that consume excluded substances
		val newDesiredRates = desiredRates filter {
			case (reaction, _) => !reaction.reactants.exists { case (f, sub) => excluded.contains(sub) }
		} map {
			case (reaction, rate) => (reaction, rate * (1-minQ))
		}
		if(newDesiredRates.isEmpty)
			return actualRates

		// recurse
		val remainingRates = _computeActualRates(remainingConc, newDesiredRates)

		// combine the results
		val it = remainingRates.iterator
		var next = it.next
		for((reaction, aRate) <- actualRates) yield {
			if(reaction == next._1) {
				val rate = aRate + next._2
				if(it.hasNext)
					next = it.next
				else
					next = (null, 0)
				(reaction, rate)
			} else
				(reaction, aRate)
		}
	}

	private def computeDesiredRates(
			concentrations: Map[Any, Double],
			reactions: Iterable[OneWayReaction]
	): Iterable[(OneWayReaction, Double)] = {
		for(reaction <- reactions)
			yield (reaction, reaction.rate * reaction.reactants.map {
					case (f, sub) => (pow(concentrations(sub), f))
				}.foldLeft(1.0)(
					(a,b) => a*b
				)
			)
	}
}

/**
 * Represents a simple (non-catalysed) reaction.
 * @param reactants list of reactants along with their multiplicity
 * @param products list of products along with their multiplicity
 * @param rate the rate of the reaction
 * @param backRate the rate of the backward reaction. Optional.
 *
 * @author Tomas Mikula
 */
class Reaction(reactants: Iterable[(Int,Any)], products: Iterable[(Int,Any)], rate: Double, backRate: Double = 0) {
	if(rate <= 0)
		throw new IllegalArgumentException("Invalid rate: " + rate + ". Must be positive.")
	if(backRate < 0)
		throw new IllegalArgumentException("Invalid backward rate: " + backRate + ". Must be non-negative.")

	private val forward = new OneWayReaction(reactants, products, rate)
	private val backward = if(backRate > 0) new OneWayReaction(products, reactants, backRate) else null

	override def toString(): String = {
		forward + ":" + backRate
	}
}

private class OneWayReaction(reacts: Iterable[(Int,Any)], prods: Iterable[(Int,Any)], r: Double) {
	def rate = r
	def reactants = reacts
	def products = prods
	def factors = (reacts map { case(f, sub)=>(-f, sub) }) ++ products

	override def toString(): String = {
		reacts.map { case (f, sub) => f + "*" + sub }.mkString("+") +
		" -> " +
		prods.map { case (f, sub) => f + "*" + sub }.mkString("+") +
		" @" + r
	}
}